### **************************************** ####
### Author: Piatinskii M., &copy 2019        ####
### License: BSD                             ####
### Source package reference:                ####
### https://github.com/cfree14/datalimited2/ ####
### **************************************** ####

# install packages if not exist yet
# install.packages("dplyr")
# install.packages("icesAdvice")
# install.packages("rmarkdown")
# install.packages("tidyverse")
# devtools::install_github("cfree14/datalimited2")

library("datalimited2")
library("dplyr")
library("tidyverse")
library("rmarkdown")
library("icesAdvice")
library("ggplot2")

# load patch from cmsy2.R
source("package/cmsy2.R")

#### ----- USER CONFIGURATION ----- ####
# define population resilience level (see fishbase or Musick, 1999)
config.population.resilience <- "Medium" # Allowed resilience level: "Very low", "Low", "Medium", "High"
config.retro.years <- 0 # retrospective analysis time window
config.population.name <- "Sprattus sprattus"
config.population.area <- "Black sea, Russian waters"
config.report.author <- "Kulba S."
config.blim_frommsy <- 0.5 # Define multiplier a to get Blim refpts, Blim = a * Bmsy
config.forecast.use = TRUE # use forecasting procedure? Configure input/forecast.csv file!
#### ----- END CONFIG SECTION ----- ####


# DO NOT TOCH ENYTHING BELOW IF YOU HAVE NO IDEA WHAT ARE YOU DOING!!!

### Some model info ###
# get population resilience information: ?resilence
# get model performance: ?performance
### /end info ###

# read model data
data <- read.csv("input/data.csv", sep = ",")
# fix cyrilic comma decimal separator
data$catch <- as.numeric(gsub(",", ".", data$catch))

if(config.forecast.use == TRUE) {
  scenarios <- read.csv("input/forecast.csv", sep = ",")
  
  # check if first year in scenarious has lag = 1 after retro data
  if (scenarios$year[1] - data$year[length(data$year)] != 1) {
    stop(paste0("Scenarious year vector has lag > 1 year after retro vector!!! First year in scenarious should be: ", data$year[length(data$year)]+1))
  }
}

# get start/terminal year from data
year.start <- data[1,1]
year.terminal <- data[length(data[,1]),1]

# fit CMSY model for actual data
cmsy <- cmsy2(year=data$year, catch=data$catch, resilience = config.population.resilience)
# fix f_hi > 1
for (i in 1:length(cmsy$ref_ts$f_hi)) {
  if (cmsy$ref_ts$f_hi[i] > 1)
    cmsy$ref_ts$f_hi[i] <- 1
}
# fix f > 1
for (i in 1:length(cmsy$ref_ts$f)) {
  if (cmsy$ref_ts$f[i] > 1)
    cmsy$ref_ts$f[i] <- 1
}
# recalc f/fmsy
cmsy$ref_ts$ffmsy <- cmsy$ref_ts$f /  cmsy$ref_pts[4, "est"]
cmsy$ref_ts$ffmsy_hi <- cmsy$ref_ts$f_hi /  cmsy$ref_pts[4, "est"]

# make forecast scenarious
if (config.forecast.use) {
  forecast <- cmsy_forecast(cmsy, scenarios)
  # build df 
  forecast.b <- data.frame(year = scenarios$year)
  forecast.bbmsy <- forecast.b
  for (s in names(forecast)) {
    sc <- forecast[[s]]
    forecast.b[[s]] <- forecast[[s]]$b
    forecast.bbmsy[[s]] <- forecast[[s]]$b / cmsy$ref_pts[5, "est"]
  }
  
  forecast.f <- data.frame(year = scenarios$year)
  forecast.ffmsy <- forecast.f
  for (s in names(forecast)) {
    sc <- forecast[[s]]
    forecast.f[[s]] <- forecast[[s]]$f
    forecast.ffmsy[[s]] <- forecast[[s]]$f / cmsy$ref_pts[4, "est"]
  }
}

# check if retro info available
if (config.retro.years == FALSE || config.retro.years < 1) {
  render("Report_noretro.Rmd", output_file = "Report.html")
} else {
  # hard fit retrospective procedure ... this can take a lot of time!
  retro.cmsy = list()
  # uhhh, lets try run cmsy2 5times ))
  for (i in 1:config.retro.years) {
    # slice data frame by -1 year
    rown <- nrow(data)+1
    data.slice <- data[(-rown+i):-rown,]
    # perform cmsy. Man, you can take a smoke at ~30 min related to your PC hardware performance ;)
    retro.cmsy[[i]] <- cmsy2(year=data.slice$year, catch=data.slice$catch, resilience = config.population.resilience)
  }
  
  # extract SSB, F from retro estimates
  retro.data <- list()
  for (i in 1:config.retro.years) {
    retro.data[[i]] <- list(
      year = retro.cmsy[[i]]$ref_ts$year,
      b = retro.cmsy[[i]]$ref_ts$b, 
      f = retro.cmsy[[i]]$ref_ts$f
    )
  }
  
  # prepare data for Mohn-rho tests
  retro.ssb.5y <- list()
  retro.f.5y <- list()
  for (i in 1:config.retro.years) {
    ssb.t <- as.data.frame(retro.data[[i]]) %>% 
      filter(year >= (year.terminal-config.retro.years)) %>% 
      select(b) %>% 
      pull(b) %>% 
      c(.,rep(NA, i))
    
    f.t <- as.data.frame(retro.data[[i]]) %>% 
      filter(year >= (year.terminal-config.retro.years)) %>% 
      select(f) %>% 
      pull(f) %>%
      c(.,rep(NA, i))
    
    retro.ssb.5y[[i]] <- ssb.t
    retro.f.5y[[i]] <- f.t
  }
  # get base year 
  retro.ssb.base <- cmsy$ref_ts %>%
    filter(year >= (year.terminal-config.retro.years)) %>%
    select(b) %>%
    pull(b)
  
  retro.f.base <- cmsy$ref_ts %>%
    filter(year >= (year.terminal-config.retro.years)) %>%
    select(f) %>%
    pull(f)
  
  # build ready-to-use mohn-rho test data.frame - B
  retro.rho.ssbdata <- list(year=(year.terminal-config.retro.years):year.terminal) %>%
    append(., list(base=retro.ssb.base)) %>%
    append(retro.ssb.5y) %>%
    as.data.frame(.)
  
  names(retro.rho.ssbdata) <- c("year", "b.base", seq(-1, -config.retro.years, by=-1))
  rownames(retro.rho.ssbdata) <- retro.rho.ssbdata[,"year"]
  retro.rho.ssbdata <- retro.rho.ssbdata[,-1]
  
  # build ready-to-use mohn-rho test data.frame - F
  retro.rho.fdata <- list(year=(year.terminal-config.retro.years):year.terminal) %>%
    append(., list(base=retro.f.base)) %>%
    append(retro.f.5y) %>%
    as.data.frame(.)
  
  names(retro.rho.fdata) <- c("year", "f.base", seq(-1, -config.retro.years, by=-1))
  rownames(retro.rho.fdata) <- retro.rho.fdata[,"year"]
  retro.rho.fdata <- retro.rho.fdata[,-1]
  
  # calc mohn-rho B
  retro.rho.ssb <- mohn(retro.rho.ssbdata, peels = config.retro.years)
  # calc mohn-rho F
  retro.rho.f <- mohn(retro.rho.fdata, peels = config.retro.years)
  
  # get max biomass value
  b.max <- max(cmsy$ref_ts$b, cmsy$ref_pts$est[3], cmsy$ref_pts$est[5]) * 1.2
  b.max.floor <- signif(b.max, digits=2)
  
  # get max F value
  f.max <- max(cmsy$ref_ts$f, cmsy$ref_pts$est[4]) * 1.2
  
  # render output report
  render("Report.Rmd")  
}
